#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <Eigen/Sparse>
#include "Poisson.hpp"
#include <json.hpp>

using json = nlohmann::json;
void ComputeConvergenceRates(const std::vector<int>& grid_resolutions, const std::vector<double>& errors) {
    std::vector<double> convergence_rates;
    for (size_t i = 1; i < grid_resolutions.size(); ++i) {
        double h1 = 1.0 / (grid_resolutions[i-1] + 1);
        double h2 = 1.0 / (grid_resolutions[i] + 1);
        double rate = std::log(errors[i-1] / errors[i]) / std::log(h1 / h2);
        convergence_rates.push_back(rate);
    }

    // Output convergence rates
    std::cout << "Convergence Rates:\n";
    for (size_t i = 0; i < convergence_rates.size(); ++i) {
        std::cout << "From " << grid_resolutions[i] << " to " << grid_resolutions[i+1] 
                  << ": " << convergence_rates[i] << "\n";
    }
}
// Function to compute the L1, L2, and L∞ norms of the error
void ComputeErrors(const std::vector<double>& U, const std::vector<double>& U_exact, double h, double& L1, double& L2, double& Linf) {
    L1 = 0.0;
    L2 = 0.0;
    Linf = 0.0;
    for (size_t i = 0; i < U.size(); ++i) {
        if (U[i] != -1) {
            double error = std::abs(U[i] - U_exact[i]);
        L1 += error;
        L2 += error * error;
        Linf = std::max(Linf, error);
        }
        
    }
    L1 *= h * h;
    L2 = std::sqrt(L2 * h * h);
}

// 将字符串转换为 DomainParams::Type
DomainParams::Type stringToDomainType(const std::string& type) {
    if (type == "SQUARE") {
        return DomainParams::SQUARE;
    } else if (type == "SQUARE_MINUS_DISK") {
        return DomainParams::SQUARE_MINUS_DISK;
    } else {
        throw std::runtime_error("Unknown domain type: " + type);
    }
}

// 将字符串转换为 BoundaryConditions::Type
BoundaryConditions::Type stringToBoundaryType(const std::string& type) {
    if (type == "DIRICHLET") {
        return BoundaryConditions::DIRICHLET;
    } else if (type == "NEUMANN") {
        return BoundaryConditions::NEUMANN;
    } else if (type == "Mixed") {
        return BoundaryConditions::Mixed;
    } else {
        throw std::runtime_error("Unknown boundary type: " + type);
    }
}


const double PI = std::acos(-1);

// 目标解 u(x, y) = sin(πx) * sin(πy)
class ExactSolution : public Function {
    public:
        double operator()(double x, double y) const override {
            return std::sin(PI*x)*std::sin(PI*y);
        }
    };
    
    class u1 : public Function {
    public:
        double operator()(double x, double y) const override {
            return 0;
        }
    };
    
    class u2 : public Function {
        public:
            double operator()(double x, double y) const override {
                return 0;
            }
        };
    
    class u3 : public Function {
        public:
            double operator()(double x, double y) const override {
                return 0;
            }
        };
    
    class u4 : public Function {
        public:
            double operator()(double x, double y) const override {
                return 0;
            }
        };
    
    // 源项 f(x, y) = -Δu = -(-sin(x) + cos^2(x) + 1) * exp(y + sin(x))
    class SourceTerm : public Function {
    public:
        double operator()(double x, double y) const override {
            double laplacian_u = -2*PI*PI*std::sin(PI*x)*std::sin(PI*y);
            return -laplacian_u;
        }
    };

class u_partial_y : public Function {
    public:
        double operator()(double x, double y) const override {
            return PI*std::cos(PI*y)*std::sin(PI*x);
        }
};

class u_partial_x : public Function {
    public:
        double operator()(double x, double y) const override {
            return PI*std::cos(PI*x)*std::sin(PI*y);
        }
};


int main() {
    // Read input file
    std::ifstream input_file("input.json");
    json input_data;
    input_file >> input_data;

    // Parse domain parameters
    DomainParams domain;
    domain.type = stringToDomainType(input_data["domain"]["type"].get<std::string>());
    domain.Point_x = input_data["domain"]["Point_x"].get<double>();
    domain.Point_y = input_data["domain"]["Point_y"].get<double>();
    domain.radius = input_data["domain"]["radius"].get<double>();

    BoundaryConditions bc;
    bc.outer_bc = stringToBoundaryType(input_data["boundary_conditions"]["outer_bc"].get<std::string>());
    bc.u_1 = std::make_shared<u1>();
    bc.u_2 = std::make_shared<u1>();
    bc.u_x = std::make_shared<u_partial_x>();  
    bc.u_y = std::make_shared<u_partial_y>(); 
    bc.inner_u_x = std::make_shared<u_partial_x>();  
    bc.inner_u_y = std::make_shared<u_partial_y>(); 
    
    bc.one_point_value = 0;
    

    // Parse source term
    auto source_term = std::make_shared<SourceTerm>();

    // Test on grids with n = 8, 16, 32, 64
    std::vector<int> grid_resolutions = input_data["grid_resolution"].get<std::vector<int>>();
    std::vector<double> L1_errors, L2_errors, Linf_errors;

    for (int n : grid_resolutions) {
        bc.outer_bc = BoundaryConditions::Mixed;
        bc.u_x = std::make_shared<u_partial_x>();  
        bc.u_y = std::make_shared<u_partial_y>();
        bc.inner_u = std::make_shared<ExactSolution>();
        bc.inner_u_x = std::make_shared<u_partial_x>();  
        bc.inner_u_y = std::make_shared<u_partial_y>();
        bc.flag_inner_u = 0;
        bc.flag_outer_u = {1, 1, 1, 1};
        bc.one_point_loc = {0,0};
        bc.one_point_value = 0.0;
        domain.grid_resolution = n;
        PoissonSolver solver(domain, bc, source_term);
        solver.Build_System();
        solver.Solve();

        const std::vector<double>& solution = solver.GetSolution();

        // Compute exact solution
        std::vector<double> U_exact((n + 2) * (n + 2), 0);
        double h = 1.0 / (n + 1);
        for (int j = 0; j < n + 2; ++j) {
            for (int i = 0; i < n + 2; ++i) {
                double x = i * h, y = j * h;
                U_exact[i + j * (n + 2)] = std::sin(PI*x)*std::sin(PI*y);
            }
        }

        // Compute errors
        double L1, L2, Linf;
        ComputeErrors(solution, U_exact, h, L1, L2, Linf);
        L1_errors.push_back(L1);
        L2_errors.push_back(L2);
        Linf_errors.push_back(Linf);

        // 生成文件名
        char filename[100];
        sprintf(filename, "Problem2_parm%d.txt", n);  // 修正文件名生成
        FILE* fp = fopen(filename, "w+");
        if (!fp) {
            std::cerr << "Failed to open file: " << filename << std::endl;
            continue;  // 如果文件打开失败，跳过当前循环
        }
        const std::vector<std::vector<int>>& map = solver.GetMap();
        // 输出数值解到文件
        for (int i = 0; i < n + 2; i++) {
            for (int j = 0; j < n + 2; j++) {
                if (map[i][j] != -1) {  // 只处理有效点
                    // std::cout << "U[" << i << "," << j << "] = " << solution[i + j * (m + 2)] << "\n";
                    fprintf(fp, "%lf %lf %lf\n", i * h, j * h, solution[i + j * (n + 2)]);
                }
            }
        }
        fclose(fp);
    }

    // Report errors and convergence rates
    std::cout << "Function2 SQUARE_MINUS_DISK NEUMANN\n";
    std::cout << "Grid Resolution\tL1 Error\tL2 Error\tL∞ Error\n";
    for (size_t i = 0; i < grid_resolutions.size(); ++i) {
        std::cout << grid_resolutions[i] << "\t" << L1_errors[i] << "\t" << L2_errors[i] << "\t" << Linf_errors[i] << "\n";
    }
    // 计算并输出收敛率
    std::cout << "L1 Error Convergence Rates:\n";
    ComputeConvergenceRates(grid_resolutions, L1_errors);

    std::cout << "\nL2 Error Convergence Rates:\n";
    ComputeConvergenceRates(grid_resolutions, L2_errors);

    std::cout << "\nL∞ Error Convergence Rates:\n";
    ComputeConvergenceRates(grid_resolutions, Linf_errors);
    return 0;
}